function [p q] = leapfrog(dt)
global U  h  g  p  q  nt dx x h1 p_global time;

dx1=1/dx;
dx2 = 1/(dx*dx);
dt2= 1/dt;

p_new=p;
q_new=q;
p_old=p;
q_old=q;
time =0;



%------ Build the matrix for the Vertical Structure variable r
b=-h*h/3*dx2;

R=zeros(2*x,2*x);
R(1,x)= b;
R(1,2)= b;
R(1,1)= -2*b;

for k=2:x-1
    R(k,k-1)=b;
    R(k,k+1) =b;
    R (k,k) = -2*b;
end

R(x,x-1)= b;
R(x,1)= b;
R(x,x)= -2*b;

r= zeros(x,1);




for n=2:nt+1; % First time step, advance using explicit euler method ?
    
    r = R_Solve(q);
    i=1;
     dpdt(i,1) = -U * (p(i+1,1) - p(x,1))/(2*dx)  - h*(q(i+1,1)-2*q(i,1)+q(x,1))/(dx*dx) + h*h/3*(r(i+1,1)-2*r(i,1)+r(x,1))/(dx*dx);
     dqdt(i,1) = -U * (q(i+1,1) - q(x,1))/(2*dx)  - g*p(i,1);
    
    for i=2:x-1    
        dpdt(i,1) = -U * (p(i+1,1) - p(i-1,1))/(2*dx)  - h*(q(i+1,1)-2*q(i,1)+q(i-1,1))/(dx*dx)+ h*h/3*(r(i+1,1)-2*r(i,1)+r(i-1,1))/(dx*dx);
        dqdt(i,1) = -U * (q(i+1,1) - q(i-1,1))/(2*dx)  - g*p(i,1);
    end
    i=x;
    dpdt(i,1) = -U * (p(1,1) - p(i-1,1))/(2*dx)  - h*(q(1,1)-2*q(i,1)+q(i-1,1))/(dx*dx)+ h*h/3*(r(1,1)-2*r(i,1)+r(i-1,1))/(dx*dx);
    dqdt(i,1) = -U * (q(1,1) - q(i-1,1))/(2*dx)  - g*p(i,1);
    
   if (n==2)
       % Use a small time step, as Explicit first order is generally
       % Unstable for Central Differenving
        p_new(1:x,1) = 0.001*dt*dpdt(1:x,1) + p(1:x,1);
        q_new(1:x,1) = 0.001*dt*dqdt(1:x,1) + q(1:x,1);

    else  
        p_new(1:x,1) = 2*dt*dpdt(1:x,1) + p_old(1:x,1);
        q_new(1:x,1) = 2*dt*dqdt(1:x,1) + q_old(1:x,1);
    end

    p_old(1:x,1) = p(1:x,1);
    q_old(1:x,1) = q(1:x,1);

    p(1:x,1)= p_new(1:x,1);
    q(1:x,1)= q_new(1:x,1);
    time =n*dt;
    if rem(time,5)==0
        k=time/5;
        p_global(:,k) = p;
        refreshdata(h1,'caller') % Evaluate p in the function workspace
        drawnow
    end

end
display('Completed Successfully');
